#include "Jacobi.h"
#include "SLAU.h"
#include "portable_time.h"
#include "function.h"
#include "stdio.h"

void main()
{

    unsigned int Sx=100, Sy=100, Sz=100;
    CRSMatrix* A=allocation(Sx,Sy,Sz);
    T dx=10.0/Sx, dy=10.0/Sy, dz=10.0/Sz, dt=1e-3; 
    T* b=new T[(Sx-2)*(Sy-2)*(Sz-2)];
    T* fu=new T[(Sx-2)*(Sy-2)*(Sz-2)];
    T* x1=new T[(Sx-2)*(Sy-2)*(Sz-2)];
    T* x2=new T[(Sx-2)*(Sy-2)*(Sz-2)];
    for (unsigned int i=0; i<Sz-2; i++)
    {
        for (unsigned int j=0; j<Sy-2; j++)
            for (unsigned int k=0; k<Sx-2; k++)
            {
                fu[k+j*(Sx-2)+i*(Sx-2)*(Sy-2)]=f(0+dx*(k+1), 0+dy*(j+1), 0+dz*(i+1));
                x1[k+j*(Sx-2)+i*(Sx-2)*(Sy-2)]=u(0+dx*(k+1), 0+dy*(j+1), 0+dz*(i+1));
            }
    }
    for (unsigned int i=0; i<(Sx-2)*(Sy-2)*(Sz-2); i++)
    {
        b[i]=fu[i]*dt+x1[i];
    }
    fillMatrix(A,1,1,1,1,dt,dx,dy,dz,Sx,Sy,Sz);
    for (unsigned int j=0; j<100; j++)
    {
        T dev=0;
        unsigned int iteration=MethodJacobi(*A,x1,b,x2,1e-9, &dev);
        T*tmp=x1;
        x1=x2;
        x2=tmp;
        for (unsigned int i=0; i<(Sx-2)*(Sy-2)*(Sz-2); i++)
        {
            b[i]=fu[i]*dt+x1[i];
        }
        printf("%i %i %e\n",j, iteration, dev);
    }
}